tsibble)와 기본 시각화(시간 플롯)
이해index(시간축)와 key(개체
식별자)를 명시한 tidy 시계열 데이터 구조이다.tsibble → feasts(진단/변환) →
fable(모형/예측)로 이어지는 tidyverts
워크플로의 출발점이다.global_economy
tourism
1998년 1분기부터 2016년 4분기까지 1박이상 호주의 인바운드 관광객 수를 나타내는 자료임
tsibble이란 객체는 R에서 다중 시계열을 저장하고 이를
제어할 수 있는 것으로 볼 수 있음
tsibble은 다음 형식을 가지고 있음:
tsibble은 tidyverse함수를 이용하여
제어한다.
tsibble
인덱스mydata <- tsibble(
year = 2015:2019,
y = c(123, 39, 78, 52, 110),
index = year
)
mydata
mydata <- tibble(
year = 2015:2019,
y = c(123, 39, 78, 52, 110)
) |>
as_tsibble(index = year)
mydata
z
z |>
mutate(Month = yearmonth(Month)) |>
as_tsibble(index = Month)
| Frequency | Function |
|---|---|
| Quarterly | yearquarter() |
| Monthly | yearmonth() |
| Weekly | yearweek() |
| Daily | as_date(), ymd() |
| Sub-daily | as_datetime() |
tibble로 변환prison <- readr::read_csv("data/prison_population.csv")
## Rows: 3072 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): state, gender, legal, indigenous
## dbl (1): count
## date (1): date
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
prison <- readr::read_csv("data/prison_population.csv") |>
mutate(Quarter = yearquarter(date))
## Rows: 3072 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): state, gender, legal, indigenous
## dbl (1): count
## date (1): date
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
prison <- readr::read_csv("data/prison_population.csv") |>
mutate(Quarter = yearquarter(date)) |>
select(-date)
## Rows: 3072 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): state, gender, legal, indigenous
## dbl (1): count
## date (1): date
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
prison <- readr::read_csv("data/prison_population.csv") |>
mutate(Quarter = yearquarter(date)) |>
select(-date) |>
as_tsibble(
index = Quarter,
key = c(state, gender, legal, indigenous)
)
## Rows: 3072 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): state, gender, legal, indigenous
## dbl (1): count
## date (1): date
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
PBS
filter()함수를 이용하여 원하는 열을 선택할 수
있다.PBS |>
filter(ATC2 == "A10")
select()함수를 이용하여 원하는 열을 선택할 수
있다.PBS |>
filter(ATC2 == "A10") |>
select(Month, Concession, Type, Cost)
summarise()함수를 이용하여 주요 변수의 요약통계량을
구할 수 있다.PBS |>
filter(ATC2 == "A10") |>
select(Month, Concession, Type, Cost) |>
summarise(TotalC = sum(Cost))
mutate() 함수를 이용해 새로운 변수를 만들어 낼 수도
있다.PBS |>
filter(ATC2 == "A10") |>
select(Month, Concession, Type, Cost) |>
summarise(TotalC = sum(Cost)) |>
mutate(Cost = TotalC / 1e6) -> a10
autoplot()/ggplot2를 활용하면 여러
시계열을 겹쳐 비교하거나, 변환(로그/차분) 전후를 손쉽게
시각화할 수 있다.a10 |>
autoplot(Cost) +
labs(y = "$ (millions)", title = "호주 항생제 판매액")
ansett |>
autoplot(Passengers)
ansett |>
filter(Class == "Economy") |>
autoplot(Passengers)
ansett |>
filter(Airports == "MEL-SYD") |>
autoplot(Passengers)
예: \(T_t\)를 추세라 하면, 대표적 표현은 선형/다항/스무딩 추세
\[ T_t = \beta_0 + \beta_1 t \quad\text{(선형)},\qquad T_t = \sum_{j=0}^p \beta_j t^j \quad\text{(다항)} \]
또는 로컬 레벨/로컬 선형(상태공간) 등 확률적 추세.
계절 주기 \(s\)에 대해 완전 주기성:
\[ S_{t+s} = S_t,\qquad \sum_{j=1}^{s} S_{t-j+1} = 0 \;(\text{식별 제약, 가산형}) \]
더미로 표현: \(S_t = \sum_{j=1}^{s} \gamma_j D_{j,t},\ \sum_{j=1}^s \gamma_j = 0\)
조화형으로도 표현 가능:
\[ S_t = \sum_{k=1}^{\lfloor s/2\rfloor} \left[a_k \cos\!\left(\tfrac{2\pi kt}{s}\right)+ b_k \sin\!\left(\tfrac{2\pi kt}{s}\right)\right] \]
감쇠 진동을 갖는 AR(2)로 근사 가능(복소근):
\[ C_t = 2\rho\cos(\omega)\,C_{t-1} - \rho^2 C_{t-2} + \varepsilon_t,\qquad \text{평균 주기} \approx \frac{2\pi}{\omega},\ \ 0<\rho<1 \]
- 시계열의 분해
가장 흔한 분해방법은 다음과 같다.
가산형(additive):
\[ y_t = T_t + S_t + C_t + E_t \]
승수형(multiplicative):
\[ y_t = T_t \times S_t \times C_t \times \varepsilon_t \]
(로그 변환 \(\tilde y_t=\log y_t\)을 쓰면 가산형으로 변환되어 추정이 용이)
여기서 \(E_t\) 또는 \(\varepsilon_t\)는 평균 0, 상수분산을 갖는 교란(통상 백색잡음 가정).
주기 길이
평균 길이
변동 폭
변환/진단
스펙트럼 관점 * 계절성은 계절 주파수 \(\tfrac{2\pi k}{s}\)에서 예리한 피크(선 스펙트럼)가 나타나는 반면, * 순환은 저주파 대역에서 완만한 능선(ridge) 형태로 에너지가 분포하는 경향.
aus_production |>
filter(year(Quarter) >= 1980) |>
autoplot(Electricity) +
labs(y = "GWh", title = "호주 전력생산량량")
aus_production |>
autoplot(Bricks) +
labs(y = "백만개개", title = "호주 벽돌 생산량")
us_employment |>
filter(Title == "Retail Trade", year(Month) >= 1980) |>
autoplot(Employed / 1e3) +
labs(y = "백만명", title = "미국 소매업 취업자수")
gafa_stock |>
filter(Symbol == "AMZN", year(Date) >= 2018) |>
autoplot(Close) +
labs(y = "미 달러", title = "아마존 주식 종가")
pelt |>
autoplot(Lynx) +
labs(y = "마리", title = "캐나다 연간 여우우 포획량")
계절성과 순환의 차이점
s=4, 월 s=12), 순환
패턴은 가변 길이.💡 핵심 계절 자료에서는 **정점(peak)**과 **저점(trough)**의 시점이 예측 가능하지만, 순환은 장기적으로 시점이 불확실함
1) 계절성의 정의와 성질
주기 \(s\)에 대해 완전 주기성:
\[ S_{t+s}=S_t,\qquad \sum_{j=1}^{s} S_{t-j+1}=0 \quad(\text{가산형 식별 제약}) \]
계절차분 \(\nabla_s=(1-B^s)\)로 제거 가능: \(\ \nabla_s y_t = y_t - y_{t-s}\).
ACF 서명: 시차 \(s, 2s, 3s,\dots\)에서 뚜렷한 스파이크.
스펙트럼: \(\tfrac{2\pi k}{s}\) (정수 \(k\))에서 선형(날카로운) 피크.
2) 순환의 정의와 근사 표현
고정 주기가 아님: 경기 국면에 따라 길이·진폭이 변함(보통 2년 이상).
**감쇠 진동 AR(2)**로 근사(복소근 존재):
\[ C_t = 2\rho\cos(\omega)\,C_{t-1} - \rho^2 C_{t-2} + \varepsilon_t,\quad 0<\rho<1, \]
\[ \text{평균 주기} \approx \frac{2\pi}{\omega},\ \ \text{진폭은 } \rho^h \text{로 감쇠} \]
ACF 서명: 사인파처럼 진동하며 서서히 감쇠(damped sinusoid).
스펙트럼: 저주파 대역에서 넓은 봉우리(능선) 형태.
3) 분해 모형과 실무 팁
4) 식별 체크리스트
계절성
순환
5) 예측 관점
계절 도표는 같은 달(또는 분기)끼리 묶어 계절 패턴을 한눈에 보여줌
계절성의 위상(peak/low), 변동성 변화(분산 상승/하락)를 쉽게 파악
a10 시계열
a10 |>
autoplot(Cost)
a10 |> gg_season(Cost, labels = "both") +
labs(y = "$ million", title = "Seasonal plot: antidiabetic drug sales")
beer <- aus_production |>
select(Quarter, Beer) |> filter(year(Quarter) >= 1992)
beer |> autoplot(Beer) +
labs(title = "호주 맥주 생산량", y = "메가리터")
beer |> autoplot(Beer) + geom_point() +
labs(title = "호주 맥주 생산량", y = "메가리터")
beer |> gg_season(Beer, labels = "right")
vic_elec
vic_elec |> autoplot()
## Plot variable not specified, automatically selected `.vars = Demand`
vic_elec |> gg_season(Demand)
vic_elec |> gg_season(Demand, period = "week")
vic_elec |> gg_season(Demand, period = "day")
a10 |>
gg_subseries(Cost) +
labs(y = "백만 달러", title = "Subseries 도표: 항생제 매출")
gg_subseries()beer <- aus_production |>
select(Quarter, Beer) |>
filter(year(Quarter) >= 1992)
beer |> autoplot(Beer)
beer |> gg_subseries(Beer)
holidays <- tourism |>
filter(Purpose == "Holiday") |>
group_by(State) |>
summarise(Trips = sum(Trips))
holidays |> autoplot(Trips) +
labs(y = "천 번", title = "호주의 국내휴가")
holidays |> gg_season(Trips) +
facet_wrap(vars(State), nrow = 2, scales = "free_y")+
labs(y = "천 번", title = "호주의 국내휴가")
holidays |>
gg_subseries(Trips) +
labs(y = "천 번", title = "호주의 국내휴가")
vic_elec_day_type <- vic_elec |>
filter(year(Time) == 2014) |>
mutate(Day_Type = case_when(
Holiday ~ "Holiday",
wday(Date) %in% 2:6 ~ "Weekday",
TRUE ~ "Weekend"))
vic_elec_day_type
vic_elec_day_type |>
ggplot(aes(x = Temperature, y = Demand)) +
geom_point() +
labs(x = "온도(섭씨)", y = "전기수요 (GW)")
vic_elec_day_type |>
ggplot(aes(x = Temperature, y = Demand, colour = Day_Type)) +
geom_point() +
labs(x = "온도(섭씨)", y = "전기수요 (GW)")
두 변수(\(y\) and \(x\)) 사이의 선형관계의 정도를 측정 .
us_change |> GGally::ggpairs(columns = 2:6)
new_production <- aus_production |>
filter(year(Quarter) >= 1992)
new_production
new_production |> gg_lag(Beer, geom = "point")
new_production |> gg_lag(Beer)
new_production |> gg_lag(Beer, geom = "point")
표본 자기공분산(lag \(k\))을 \(c_k\), 표본 자기상관을 \(r_k\)라 하자. 표본평균 \(\bar y=\frac1T\sum_{t=1}^T y_t\)일 때,
\[ c_k \;=\; \frac{1}{T}\sum_{t=k+1}^T (y_t-\bar{y})(y_{t-k}-\bar{y}), \qquad r_k \;=\; \frac{c_k}{c_0}. \]
참고: 일부 교과서/소프트웨어는 \(c_k\) 분모를 \(T-k\)로 두는 불편(unbiased) 정의를 쓰기도 한다(위 식은 bias 버전). 핵심 개념과 해석은 동일하다.
첫 9개 시차의 자기상관을 계산:
new_production |> ACF(Beer, lag_max = 9)
시각화(상관도표, correlogram):
new_production |> ACF(Beer, lag_max = 9) |> autoplot()
용어
여러 시차의 \(r_k\)들을 모은 벡터 \((r_1,r_2,\dots)\)를 ACF(autocorrelation function)라 한다.
이를 막대그래프로 그린 것을 상관도표(correlogram) 라고 한다.
new_production |> ACF(Beer) |> autoplot()
ACF: 전체 자기상관(중간 시차 효과 포함)
PACF: 부분 자기상관(중간 시차를 통제한 순수 효과)
1) 왜 추세가 ACF를 크게 만들까?
2) 추세 제거(차분) 후 ACF 해석
3) 계절성의 ACF 서명과 제거
4) 진단·실무 팁
# 추세·계절성 진단과 차분 예시 (feasts/tidyverts)
library(tsibble)
library(feasts)
library(dplyr)
# 예: tsibble 객체 'x'와 변수 'y'가 있다고 가정
x |> ACF(y) |> autoplot() # 원자료 ACF
# 1차 차분 후
x |> mutate(dy = difference(y)) |>
ACF(dy) |> autoplot()
# 계절 차분(s=4) + 1차 차분 후
x |> mutate(dsy = difference(y, lag = 4),
dsy1 = difference(dsy)) |>
ACF(dsy1) |> autoplot()
retail <- us_employment |>
filter(Title == "Retail Trade", year(Month) >= 1980)
retail |> autoplot(Employed)
retail |>
ACF(Employed, lag_max = 48) |>
autoplot()
google_2015 <- gafa_stock |>
filter(Symbol == "GOOG", year(Date) == 2015) |>
select(Date, Close)
google_2015
google_2015 |> autoplot(Close)
google_2015 |>
ACF(Close, lag_max = 100)
google_2015 |>
ACF(Close, lag_max = 100) |>
autoplot()
set.seed(30)
wn <- tsibble(t = 1:50, y = rnorm(50), index = t)
autoplot(wn, y) + labs(title = "백색잡음 시계열 예시", x = "t", y = "y_t")
정의(이론)
백색잡음 \(\{y_t\}\)는 통상 다음을 가정함 \[ \mathbb{E}(y_t)=0,\qquad \mathrm{Var}(y_t)=\sigma^2\ \text{(상수)},\qquad \mathrm{Cov}(y_t,y_{t-k})=0\ \ (k\ne 0). \] 실무에서는 “시간에 따라 상관이 없고 평균 0, 분산이 일정”하면 백색잡음으로 보기도 함 (엄밀히는 독립성까지 요구하는 경우가 많음.)
| \(r_{1}\) | \(r_{2}\) | \(r_{3}\) | \(r_{4}\) | \(r_{5}\) | \(r_{6}\) | \(r_{7}\) | \(r_{8}\) | \(r_{9}\) | \(r_{10}\) |
|---|---|---|---|---|---|---|---|---|---|
| 0.014 | -0.163 | 0.163 | -0.259 | -0.198 | 0.064 | -0.139 | -0.032 | 0.199 | -0.024 |
wn |>
ACF(y) |>
autoplot() +
labs(title = "백색잡음의 ACF", x = "Lag", y = "ACF")
표본분포(이론 보충)
백색잡음에서 \(r_k\)는 큰 표본에서 \[ r_k \ \overset{approx}{\sim}\ \mathcal{N}\!\left(0,\ \frac{1}{T}\right), \] 이므로 약 95%의 \(r_k\)가 \(\pm 1.96/\sqrt{T}\) 안에 들 것으로 기대함
여러 시차를 동시에 검정할 때는 Ljung–Box 검정을 활용
# Ljung-Box: 잔차의 백색잡음 여부 테스트(여기서는 자체 시계열)
Box.test(wn$y, lag = 10, type = "Ljung-Box", fitdf = 0)
##
## Box-Ljung test
##
## data: wn$y
## X-squared = 12.993, df = 10, p-value = 0.2241
pigs <- aus_livestock |>
filter(State == "Victoria", Animal == "Pigs", year(Month) >= 2014)
autoplot(pigs, Count/1e3) +
labs(y = "천 마리(Thousands)",
title = "빅토리아주 월별 돼지 도축 수 (2014–2018)")
pigs |>
ACF(Count) |>
autoplot() +
labs(title = "돼지 도축수의 ACF", x = "Lag(월)", y = "ACF")
모형화 힌트(이론 보충)
- 계절성이 보이면 계절 차분 \(\nabla_{12} y_t = y_t - y_{t-12}\) 또는
계절 모형화(SARIMA, ETS 계절 성분)를 고려
- 추세와 계절이 함께 있을 때는 \(\nabla^d
\nabla_{12}^D y_t\)처럼 복합 차분을 적용해
정상화 후 ACF/PACF로 차수를 판단
- 백색잡음/잔차 진단에는 ACF/PACF와 함께 Ljung–Box 검정을 병행함